还是用RSeQC对比对后的转录组数据做一下质控

早在3年前我就关注过这个软件，因为某个公司的转录组结题报告里面提到过，但那个时候我还是不是很懂NGS组学异同，不太理解质量控制，只知道流程，别人分析什么我就分析什么，所以那个时候我博客里面关于RSeQC的教程很有趣(http://www.bio-info-trainee.com/113.html) 复制粘贴该url到浏览器可以查看我当初的教程，或者点击最下面的阅读原文即可。

那个时候写教程，以软件安装，软件input和output为主，因为觉得新手最容易纠结的就是这些了，但是现在回过头来看，软件安装已经成了小菜一碟，对各种bam/sam/vcf/gtf也耳熟能详，所谓的input/output也不是问题了。

所以，再看看我最近是如何记录该软件的吧：

RSeQC包是一个python软件，最新版是 v2.6.4 ， 依赖于： gcc; python2.7; numpy; R

它提供了一系列有用的小工具能够评估高通量测序尤其是RNA-seq数据，比如一些基本模块，检查序列质量, 核酸组分偏性, PCR偏性, GC含量偏性,还有RNA-seq特异性模块: 评估测序饱和度， 映射读数分布， 覆盖均匀性， 链特异性， 转录水平RNA完整性等

详细列表如下：

• bam2fq.py

• bam2wig.py

• bam_stat.py

• clipping_profile.py

• deletion_profile.py

• divide_bam.py

• FPKM_count.py

• geneBody_coverage.py

• geneBody_coverage2.py

• infer_experiment.py

• inner_distance.py

• insertion_profile.py

• junction_annotation.py

• junction_saturation.py

• mismatchprofile.pynormalizebigwig.pyoverlay_bigwig.py

• read_distribution.py

• read_duplication.py

• readGC.pyreadhexamer.py

• read_NVC.py

• read_quality.py

• RNAfragmentsize.py

• RPKMcount.pyRPKMsaturation.py

• spilt_bam.py

• splitpairedbam.py

• tin.py

数据库文件

RSeQC接受4种文件格式:

• BED 格式: Tab 分割, 12列的表示基因模型的纯文本文件

• SAM 或BAM 格式: 用来存储reads 比对结果信息.

• 染色体大小文件: 只有两列的纯文本文

• Fasta文件的参考基因组

数据库文件根据参考基因组版本自行选择下载，我这里要下载的是hg19系列，下载地址如下：

1. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_GENCODE_v19_basic.Cancer_genes.bed.gz

2. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_UCSC_knownGene.bed.gz

3. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_RefSeq.bed.gz

4. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_rRNA.bed.gz

5. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_GENCODE_GENE_V19_comprehensive.bed.gz

6. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19.HouseKeepingGenes.bed

7. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_AceView.bed.gz

8. https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_Ensembl.bed.gz

9. ls *.gz|while read id ; do gunzip $id;done

希望读者能够明白，看教程一定要看规律，我为什么列出如此多的url，其实就是想你领悟它们的共性： https://sourceforge.net/projects/rseqc/files/ 你在浏览器打开就明白了。

### 软件安装

1. # 如果python版本没有问题，那么直接用pip即可安装

2. pip install RSeQC --user

3. # 如果有conda，那么更方便

4. conda install -c bioconda rseqc

5. ## 依赖于python2.7

6. ## 所以conda可能需要先创建python2.7的环境,再安装

7. conda info --envs

8. conda create -n py2.7 python=2.7 rseqc

9. source activate py2.7

虽然该软件的使用命令非常多，但很多功能并不是用来诊断转录组测序的，所以不在我们的考虑范围内。下面是我们经常会用得到的：

1. # nohup  bash run.sh 1>run.log 2>&1 &

2. #source activate py2.7  

3. mkdir -p db

4. cd db

5. if [ ! -f hg19_RefSeq.bed ]; then

6. wget https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_RefSeq.bed.gz

7. wget https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg38_RefSeq.bed.gz

8. wget https://sourceforge.net/projects/rseqc/files/BED/Mouse_Mus_musculus/mm10_RefSeq.bed.gz

9. ls *.gz|while read id ; do gunzip $id;done

10. fi

11. cd ../

12. bed='db/mm10_RefSeq.bed'

13. nohup geneBody_coverage.py -r $bed -i bam_path.txt  -o coverage 1>coverage.log 2>&1 &

14. cat bam_path.txt |while read id ; do

15.      echo $id

16.      file=$(basename $id )

17.      sample=${file%%.*}

18.      junction_annotation.py -i   $id   -r $bed  -o  ${sample}_junction

19.      bam_stat.py -i $id >${sample}_bam_stat.log  

20.      RPKM_saturation.py -r $bed -d '1++,1--,2+-,2-+' -i $id -o ${sample}_RPKM_saturation

21.      ## below just like fastqc

22.      nohup read_quality.py -i  $id -o ${sample}_read_quality &

23.      nohup read_NVC.py -i  $id -o ${sample}_read_NVC &

24.      nohup read_GC.py -i $id -o ${sample}_read_GC &

25.      nohup read_duplication.py -i $id -o ${sample}_read_duplication &

26.      read_distribution.py -i  $id  -r $bed  >${sample}_distribution.log

27. done  

用 bam_stat.py来统计总比对记录, PCR重复数, Non Primary Hits表示多匹配位点, 不匹配的reads数, 比对到+链的reads, 比对到-链的reads, 有剪切位点的reads等.

1. #==================================================

2. #All numbers are READ count

3. #==================================================

4. Total records:                          45722327

5. QC failed:                              0

6. Optical/PCR duplicate:                  0

7. Non primary hits                        2687735

8. Unmapped reads:                         2338796

9. mapq < mapq_cut (non-unique):           2045264

10. mapq >= mapq_cut (unique):              38650532

11. Read-1:                                 19631272

12. Read-2:                                 19019260

13. Reads map to '+':                       19320271

14. Reads map to '-':                       19330261

15. Non-splice reads:                       20690614

16. Splice reads:                           17959918

17. Reads mapped in proper pairs:           36737552

18. Proper-paired reads map to different chrom:0

可以看到比对效果非常赞，这个转录组很成功！

另外一个比较赞的小程序就是： read_duplication.py 结果一般如下：

1. Total Reads                   40695796

2. Total Tags                    64718115

3. Total Assigned Tags           61411678

4. =====================================================================

5. Group               Total_bases         Tag_count           Tags/Kb

6. CDS_Exons           34406515            45257520            1315.38

7. 5'UTR_Exons         6859302             2274659             331.62

8. 3'UTR_Exons         25952114            9778098             376.77

9. Introns             943281009           3254031             3.45

10. TSS_up_1kb          19391072            65573               3.38

11. TSS_up_5kb          88202906            155561              1.76

12. TSS_up_10kb         160360035           222457              1.39

13. TES_down_1kb        19659116            216878              11.03

14. TES_down_5kb        84349049            524626              6.22

15. TES_down_10kb       149723035           624913              4.17

16. =====================================================================

可以用一个饼图来表示，在生信技能树论坛里面还有人专门提问过。

用 geneBody_coverage.py来计算RNA-seq reads在基因上的覆盖度，这里推荐对所有的样本的 bam文件一起运行该程序进行诊断，如图：

junction_annotation.py:

输入一个 BAM或 SAM文件和一个 bed12格式的参考基因文件,这个模块将根据参考基因模型计算剪切融合(splice junctions)事件.

• splice read: 一个RNA read,能够被剪切一次或多次

• splice junction:多个跨越同一个内含子的剪切事件能够合并为一个 splicing junction.

一般来说，novel的junction区域总是有的，因为我们用的是ucsc的refseq参考注释集，本身就是不够完整的。

RPKM_saturation.py

  任何样本统计（ RPKM）的精度受样本大小（ 测序深度）的影响，重抽样或切片是使用部分数据来评估样本统计量的精度的方法。这个模块从总的 RNA reads中重抽样并计算每次的 RPKM值，通过这样我们就能检测当前测序深度是不是够的(如果测序深度不够RPKM的值将不稳定,如果测序深度足够则RPKM值将稳定)。*默认情况下,这个模块将计算20个 RPKM值(分别是对个转录本使用5%,10%,…,95%的总 reads)，所以非常消耗内存哦。

  在结果图中,Y轴表示 “Percent Relative Error” 或 “Percent Error”

说明:Q1,Q2,Q3,Q4是按照转录本表达量4分位分开的.Q1表示的是表达量低于25%的转录本,以此类推.可以看出:随着样本量升高, RPKM与实际值的偏差也在降低.而且转录本表达量越高这种趋势越明显(Q4最明显).

写在最后：

NGS组学分析流程的每一个步骤都应该是有充分的质量控制，主要是考虑到各个项目的实际情况可能会比较特殊，如果都走一样的自动化，流水线的流程，肯定是会有问题的。

明天给大家看看，问题主要是什么，敬请期待哈。